GLMM example 4

Preparations

Load the necessary libraries

library(car) # for regression diagnostics
library(broom) # for tidy output
library(broom.mixed)
library(ggfortify) # for model diagnostics
library(sjPlot) # for outputs
library(knitr) # for kable
library(effects) # for partial effects plots
library(ggeffects) # for effects plots in ggplot
library(emmeans) # for estimating marginal means
library(MASS) # for glm.nb
library(MuMIn) # for AICc
library(tidyverse) # for data wrangling
library(DHARMa) # for residuals and diagnostics
library(nlme) # for lme
library(lme4) # for glmer
library(glmmTMB) # for glmmTMB
library(performance) # for diagnostic plots
library(see) # for diagnostic plots

Scenario

Crab_shrimp_coral

To investigate synergistic coral defence by mutualist crustaceans, Mckeon et al. (2012) conducted an aquaria experiment in which colonies of a coral species were placed in a tank along with a predatory sea star and one of four symbiont combinations:

• no symbiont,
• a crab symbiont
• a shrimp symbiont
• both a crab and shrimp symbiont.

The experiments were conducted in a large octagonal flow-through seawater tank that was partitioned into eight sections, which thereby permitted two of each of the four symbiont combinations to be observed concurrently. The tank was left overnight and in the morning, the presence of feeding scars on each coral colony was scored as evidence of predation. The experiments were repeated ten times, each time with fresh coral colonies, sea stars and symbiont.

The ten experimental times represent blocks (random effects) within which the symbiont type (fixed effect) are nested.

Read in the data

mckeon <- read_csv("../public/data/mckeon.csv", trim_ws = TRUE)
glimpse(mckeon)

## Rows: 80
## Columns: 3
## $ BLOCK     <dbl> 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10…
## $ PREDATION <dbl> 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, …
## $ SYMBIONT  <chr> "none", "none", "none", "none", "none", "none", "none", "non…

Since the response here is the presence or absence of predation (feeding scars), a binomial distribution is appropriate.

We need to make sure that the categorical variable and the random effect are defined as factors. When doing so, it might be valuable to rearrange the order of the fixed effect (SYMBIONT) such that the none group is considered the first group. This way, the other levels will all naturally be compared to this level (hence it will be treated as a reference of control group).

mckeon <- mckeon %>%
  mutate(
    BLOCK = factor(BLOCK),
    SYMBIONT = factor(SYMBIONT, levels = c("none", "crabs", "shrimp", "both"))
  )

Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{N}(n, p_i)\\ ln\left(\frac{p_i}{1-p_1}\right) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]

where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of symbionts on the probability of the colony experiencing predation. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual coral colonies.

ggplot(mckeon, aes(y = PREDATION, x = SYMBIONT)) +
  geom_point(position = position_jitter(width = 0.2, height = 0)) +
  facet_wrap(~BLOCK)

Fit the model

The data have been setup as a single factor (mixed effects) design. Alternatively, we could have coded this up as a factorial (mixed effects) such that there was a dummy variable for Crabs, and a dummy variable for Shrimps. We could then fit a model that had the main effects of Crabs and Shrimp as well as there interaction.

Since models with categorical variables essentially dummy code up the variable anyway, fitting a single factor model will result in similar parameter estimates to a factorial design with the exception of the interaction (or Crabs and Shrimp group). Nevertheless, having fit the model as a single factor model, we can then perform a specific contrast that is the same as the interaction effect from a factorial design.

glmer (lme4)

Lets start with the random intercept model

mckeon.glmer1 <- glmer(PREDATION ~ SYMBIONT + (1 | BLOCK),
  data = mckeon,
  family = binomial(link = "logit")
)

And now attempt the random intercept/slope model. By default, glmer uses the bobyqa optimiser. This optimiser uses iterative quadratic approximation of twice differentiable functions. When functions are not twice differentiable, it may perform poorly.

We could specify the new random intercept/slope model out in full or we could simply use the update() function to make a modification (in this case, we can remove the random intercept and add the random intercept/slope).

mckeon.glmer2 <- update(mckeon.glmer1, ~ . - (1 | BLOCK) + (SYMBIONT | BLOCK))

Along with attempting to fit the model, a glmer() performs a set of model checks to help diagnose whether the optimisation process has successfully converged as well as passed other criterion.

Evidently, the above (random intercepts/slope) model failed to converge after 10000 evaluations of the likelihood function. That is, it failed to reach a consensus on the parameter estimates. Note, this is a warning not an error. It does not necessarily mean that the model fit is not correct. Nevertheless, without further exploration of what might be causing the lack of convergence, such a model should not be used as the parameter estimates (and/or their standard errors) may be unreliable.

When receiving a convergence warning, the following potential remedies are recommended:

• alternative optimizers (such as Nelder-Mead, nlminbwrap or L-BFGS-B)

• increase the number of iterations in an attempt to provide more scope for consensus

• abandon the model

• use allFit() to attempt to fit will all the available optimizers. Different optimizers perform differently under different conditions. The default optimiser (bobyqa) has been selected by the lme4 team as they consider it to be the best all-round choice and a good compromise between robustness and efficiency. However, other optimizers might work better for any given circumstance. Other common optimizers include:

  • Nelder-Mead: This is a downhill Simplex method
  • non-linear minimisation with box constraint (nlminb)
  • large-scale, quasi-Newton bound-constrained (L-BCGS-B) method

  If all of the optimizers yield very similar parameter estimates, then the convergence warnings can be considered as false positives and as such, ignored (unless they all suffer the same structural issue due to a poorly specified or over-fit model).

• adjust the stopping tolerance for the optimiser (e.g. optCtrl=list(maxfun=2e5))

• center (and possibly scale) continuous predictors

• consider a more expensive Hessian calculation. Hessian computations are used for convergence checking as well as estimating the standard errors of fixed effects parameters in GLMM. By default, the method used provides very rapid estimates. However, they can be unreliable.

• restart the fit with starting values taken from the current parameter estimates

allFit

Let start with the first of these options - often it is the best initial approach.

mckeon.allFit <- mckeon.glmer2 %>% allFit()

## bobyqa : [OK]
## Nelder_Mead : [OK]
## nlminbwrap : [OK]
## nmkbw : [OK]
## optimx.L-BFGS-B : [OK]
## nloptwrap.NLOPT_LN_NELDERMEAD : [OK]
## nloptwrap.NLOPT_LN_BOBYQA : [OK]

## Check which of the models are considered valid (OK)
is.OK <- sapply(mckeon.allFit, is, "merMod")
is.OK

##                        bobyqa                   Nelder_Mead 
##                          TRUE                          TRUE 
##                    nlminbwrap                         nmkbw 
##                          TRUE                         FALSE 
##               optimx.L-BFGS-B nloptwrap.NLOPT_LN_NELDERMEAD 
##                          TRUE                          TRUE 
##     nloptwrap.NLOPT_LN_BOBYQA 
##                          TRUE

Conclusions:

• all but the nmkbw model are considered valid (OK)

We will check the convergence characteristics of only the valid models

diff_optims.OK <- mckeon.allFit[is.OK]
lapply(diff_optims.OK, function(x) x@optinfo$conv$lme4$messages)

## $bobyqa
## [1] "boundary (singular) fit: see ?isSingular"
## 
## $Nelder_Mead
## [1] "boundary (singular) fit: see ?isSingular"
## 
## $nlminbwrap
## [1] "boundary (singular) fit: see ?isSingular"
## 
## $`optimx.L-BFGS-B`
## [1] "boundary (singular) fit: see ?isSingular"
## 
## $nloptwrap.NLOPT_LN_NELDERMEAD
## [1] "unable to evaluate scaled gradient"                                       
## [2] "Model failed to converge: degenerate  Hessian with 1 negative eigenvalues"
## 
## $nloptwrap.NLOPT_LN_BOBYQA
## [1] "boundary (singular) fit: see ?isSingular"

Conclusions:

• each of the alternative optimisation routines yielded a similar outcome - singularity.

Singularity refers to presence of one or more parameters whose estimated values are on the boundary of the possible parameter space (range of possible values). Typically, this is either that there are either variances that are estimated to be 0, or correlations that are estimated to be either -1 or 1.

Singular fits usually suggest over-fitting and are more likely to contain numerical problems and lead to issues concerning testing hypotheses on boundaries.

Situations that are likely to contribute to singularity include:

• small (<5) number of random effects levels
• complex random effects structures (such as f|g - particularly where f is a categorical variable with a relatively large number of levels)

Advice on dealing with singularity issues varies widely from simplifying the random effects structure to fitting the model in a Bayesian framework in which singularity can be avoided via appropriate priors.

For the current example, we might consider dropping the random intercept/slope.

Before excepting that we must drop the random intercept/slope, it would be worth quickly exploring the other options listed above.

Number of evaluations

Lets see if extending the number of evaluations will result in convergence.

mckeon.glmer2 <- mckeon.glmer1 %>% update(~ . - (1 | BLOCK) + (SYMBIONT | BLOCK),
  control = glmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 2e5)
  )
)

Conclusions:

• apparently not..

Center/scale predictors

In the current example, the fixed effect (SYMBIONT) is a categorical variable. Consequently, it is effectively already scaled, so the centering (and/or scaling) option isn’t a solution here.

Hessian computation

Now lets consider the more expensive Hessian computation.

pars <- getME(mckeon.glmer2, c("theta", "fixef"))
devfun <- update(mckeon.glmer2, devFunOnly = TRUE)
library(numDeriv)
cat("hess:\n")

## hess:

print(hess <- hessian(devfun, unlist(pars)))

##                [,1]          [,2]          [,3]          [,4]          [,5]
##  [1,]  0.4213975722 -2.026547e-04 -1.644069e-01  1.526993e-01 -2.826724e-04
##  [2,] -0.0002026547  1.473090e-05 -1.103229e-04  2.360222e-07  5.889714e-05
##  [3,] -0.1644068830 -1.103229e-04  5.610199e-01 -4.732408e-01 -3.033088e-04
##  [4,]  0.1526992889  2.360222e-07 -4.732408e-01  4.070909e-01 -7.901643e-05
##  [5,] -0.0002826724  5.889714e-05 -3.033088e-04 -7.901643e-05  2.743109e-04
##  [6,]  0.0595859600 -3.033426e-04 -1.002367e-01  1.109353e-01 -1.446770e-03
##  [7,] -0.0128592645 -7.902433e-05  1.109353e-01 -8.059782e-02 -3.480331e-04
##  [8,] -0.0001698081 -9.029730e-09 -1.602309e-07 -3.028307e-07 -3.996483e-09
##  [9,] -0.0010049319 -2.569394e-08 -6.284820e-07 -4.679696e-07 -1.400572e-08
## [10,] -0.0009648772 -5.504486e-08 -6.568830e-07 -3.542815e-07 -1.528416e-08
## [11,]  0.0017816129  8.246683e-04 -1.147059e-02 -2.284188e-03  3.813384e-03
## [12,]  0.0005300320 -9.729221e-05  5.407133e-04  1.195055e-04 -4.498222e-04
## [13,] -0.0353598093  7.565127e-04 -3.488098e-02 -8.343249e-03  3.498610e-03
## [14,]  0.0366325791  1.654481e-04  2.286968e-02  5.939556e-03  7.655850e-04
##                [,6]          [,7]          [,8]          [,9]         [,10]
##  [1,]  5.958596e-02 -1.285926e-02 -1.698081e-04 -1.004932e-03 -9.648772e-04
##  [2,] -3.033426e-04 -7.902433e-05 -9.029730e-09 -2.569394e-08 -5.504486e-08
##  [3,] -1.002367e-01  1.109353e-01 -1.602309e-07 -6.284820e-07 -6.568830e-07
##  [4,]  1.109353e-01 -8.059782e-02 -3.028307e-07 -4.679696e-07 -3.542815e-07
##  [5,] -1.446770e-03 -3.480331e-04 -3.996483e-09 -1.400572e-08 -1.528416e-08
##  [6,]  1.192637e-01  1.566713e-02 -1.551708e-07 -1.032991e-07 -1.750863e-07
##  [7,]  1.566713e-02  5.188491e-02  9.812057e-08 -2.048719e-07 -1.002466e-07
##  [8,] -1.551708e-07  9.812057e-08  5.836355e-01 -4.984751e-01 -1.088363e-03
##  [9,] -1.032991e-07 -2.048719e-07 -4.984751e-01  4.260798e-01 -9.365976e-04
## [10,] -1.750863e-07 -1.002466e-07 -1.088363e-03 -9.365976e-04  4.259329e-01
## [11,] -5.303501e-02 -1.056178e-02 -9.029329e-08 -4.217722e-07 -4.492473e-07
## [12,]  2.500050e-03  5.533505e-04 -1.728246e-08 -2.747050e-08 -1.480038e-08
## [13,] -1.612714e-01 -3.857528e-02 -2.337290e-07 -6.183609e-07 -5.379129e-07
## [14,]  1.057451e-01  2.746087e-02 -7.623707e-09 -8.387474e-08 -1.079613e-07
##               [,11]         [,12]         [,13]         [,14]
##  [1,]  1.781613e-03  5.300320e-04 -3.535981e-02  3.663258e-02
##  [2,]  8.246683e-04 -9.729221e-05  7.565127e-04  1.654481e-04
##  [3,] -1.147059e-02  5.407133e-04 -3.488098e-02  2.286968e-02
##  [4,] -2.284188e-03  1.195055e-04 -8.343249e-03  5.939556e-03
##  [5,]  3.813384e-03 -4.498222e-04  3.498610e-03  7.655850e-04
##  [6,] -5.303501e-02  2.500050e-03 -1.612714e-01  1.057451e-01
##  [7,] -1.056178e-02  5.533505e-04 -3.857528e-02  2.746087e-02
##  [8,] -9.029329e-08 -1.728246e-08 -2.337290e-07 -7.623707e-09
##  [9,] -4.217722e-07 -2.747050e-08 -6.183609e-07 -8.387474e-08
## [10,] -4.492473e-07 -1.480038e-08 -5.379129e-07 -1.079613e-07
## [11,]  3.960248e+00 -6.832716e-03  3.526121e-01  1.446899e-02
## [12,] -6.832716e-03  7.879640e-04 -6.203661e-03 -1.417017e-03
## [13,]  3.526121e-01 -6.203661e-03  1.164430e+00 -8.056139e-01
## [14,]  1.446899e-02 -1.417017e-03 -8.056139e-01  8.214999e-01

cat("grad:\n")

## grad:

print(grad <- grad(devfun, unlist(pars)))

##  [1]  3.633557e-03 -7.390718e-09  3.913262e-08 -5.231554e-08 -4.078292e-08
##  [6]  4.978622e-09  1.356960e-09 -6.296769e-08  5.468531e-08  1.552051e-08
## [11]  1.583982e-07 -2.453373e-08 -2.122501e-08 -5.751713e-08

cat("scaled gradient:\n")

## scaled gradient:

# print(scgrad <- solve(chol(hess), grad))
## compare with internal calculations:
mckeon.glmer2@optinfo$derivs

## $gradient
##  [1]  3.633557e-03 -7.105427e-09  4.288125e-08 -5.375256e-08 -4.074963e-08
##  [6]  4.476419e-09  1.350031e-09 -6.313172e-08  5.453415e-08  1.463718e-08
## [11]  1.515588e-07 -2.437162e-08 -2.149392e-08 -5.524470e-08
## 
## $Hessian
##                [,1]          [,2]          [,3]          [,4]          [,5]
##  [1,]  4.206285e-01 -2.052784e-04 -1.644053e-01  1.526983e-01 -2.810955e-04
##  [2,] -2.052784e-04  2.288818e-05 -1.113415e-04  7.152557e-07  5.769730e-05
##  [3,] -1.644053e-01 -1.113415e-04  5.610199e-01 -4.732451e-01 -3.056526e-04
##  [4,]  1.526983e-01  7.152557e-07 -4.732451e-01  4.070864e-01 -7.963181e-05
##  [5,] -2.810955e-04  5.769730e-05 -3.056526e-04 -7.963181e-05  2.880096e-04
##  [6,]  5.958509e-02 -3.020763e-04 -1.002369e-01  1.109383e-01 -1.446962e-03
##  [7,] -1.285934e-02 -7.915497e-05  1.109328e-01 -8.060026e-02 -3.476143e-04
##  [8,]  2.384186e-07  1.192093e-06  2.384186e-07  2.384186e-07 -7.152557e-07
##  [9,]  1.668930e-06 -1.907349e-06 -4.768372e-07 -2.384186e-07  9.536743e-07
## [10,]  1.430511e-06  2.384186e-07  1.907349e-06  1.907349e-06  9.536743e-07
## [11,]  1.781225e-03  8.265972e-04 -1.146817e-02 -2.286196e-03  3.813505e-03
## [12,]  5.316734e-04 -9.822845e-05  5.424023e-04  1.187325e-04 -4.477501e-04
## [13,] -3.535891e-02  7.569790e-04 -3.487992e-02 -8.344889e-03  3.500223e-03
## [14,]  3.663301e-02  1.678467e-04  2.287388e-02  5.942345e-03  7.662773e-04
##                [,6]          [,7]          [,8]          [,9]         [,10]
##  [1,]  5.958509e-02 -1.285934e-02  2.384186e-07  1.668930e-06  1.430511e-06
##  [2,] -3.020763e-04 -7.915497e-05  1.192093e-06 -1.907349e-06  2.384186e-07
##  [3,] -1.002369e-01  1.109328e-01  2.384186e-07 -4.768372e-07  1.907349e-06
##  [4,]  1.109383e-01 -8.060026e-02  2.384186e-07 -2.384186e-07  1.907349e-06
##  [5,] -1.446962e-03 -3.476143e-04 -7.152557e-07  9.536743e-07  9.536743e-07
##  [6,]  1.192703e-01  1.566839e-02  9.536743e-07  0.000000e+00 -7.152557e-07
##  [7,]  1.566839e-02  5.189705e-02  1.192093e-06 -1.907349e-06  2.384186e-07
##  [8,]  9.536743e-07  1.192093e-06  5.827065e-01 -4.972365e-01 -2.384186e-07
##  [9,]  0.000000e+00 -1.907349e-06 -4.972365e-01  4.245358e-01 -1.668930e-06
## [10,] -7.152557e-07  2.384186e-07 -2.384186e-07 -1.668930e-06  4.245272e-01
## [11,] -5.303502e-02 -1.056194e-02 -4.768372e-07 -4.768372e-07 -4.768372e-07
## [12,]  2.501011e-03  5.519390e-04 -2.384186e-06  2.622604e-06  7.152557e-07
## [13,] -1.612737e-01 -3.857756e-02 -2.384186e-07  2.384186e-06  1.192093e-06
## [14,]  1.057370e-01  2.746010e-02 -1.907349e-06  2.145767e-06 -4.768372e-07
##               [,11]         [,12]         [,13]         [,14]
##  [1,]  1.781225e-03  5.316734e-04 -3.535891e-02  3.663301e-02
##  [2,]  8.265972e-04 -9.822845e-05  7.569790e-04  1.678467e-04
##  [3,] -1.146817e-02  5.424023e-04 -3.487992e-02  2.287388e-02
##  [4,] -2.286196e-03  1.187325e-04 -8.344889e-03  5.942345e-03
##  [5,]  3.813505e-03 -4.477501e-04  3.500223e-03  7.662773e-04
##  [6,] -5.303502e-02  2.501011e-03 -1.612737e-01  1.057370e-01
##  [7,] -1.056194e-02  5.519390e-04 -3.857756e-02  2.746010e-02
##  [8,] -4.768372e-07 -2.384186e-06 -2.384186e-07 -1.907349e-06
##  [9,] -4.768372e-07  2.622604e-06  2.384186e-06  2.145767e-06
## [10,] -4.768372e-07  7.152557e-07  1.192093e-06 -4.768372e-07
## [11,]  3.960271e+00 -6.834507e-03  3.526089e-01  1.446557e-02
## [12,] -6.834507e-03  7.982254e-04 -6.203413e-03 -1.418829e-03
## [13,]  3.526089e-01 -6.203413e-03  1.164444e+00 -8.056104e-01
## [14,]  1.446557e-02 -1.418829e-03 -8.056104e-01  8.215151e-01

Conclusions:

• most of the gradient and Hessian values are very similar suggesting that the more complex Hessian computations would not help in this case.

Starting values

Finally, lets consider re-fitting the model using the current parameter estimates as starting values.

# strict_tol <- glmerControl(optCtrl=list(xtol_abs=1e-8, ftol_abs=1e-8))
pars <- getME(mckeon.glmer2, c("theta", "fixef"))
mckeon.glmer3 <- update(mckeon.glmer2, start = pars, control = glmerControl(
  optimizer = "bobyqa",
  optCtrl = list(maxfun = 2e5)
))

Conclusions:

• model is still singular…

glmmTBM (glmmTMB)

Lets start with the random intercept model

mckeon.glmmTMB1 <- glmmTMB(PREDATION ~ SYMBIONT + (1 | BLOCK),
  data = mckeon,
  family = binomial(link = "logit"),
  REML = TRUE
)

And now attempt the random intercept/slope model.

mckeon.glmmTMB2 <- mckeon.glmmTMB1 %>% update(~ . - (1 | BLOCK) + (SYMBIONT | BLOCK))

## Error in (function (start, objective, gradient = NULL, hessian = NULL, : NA/NaN gradient evaluation

mckeon.glmmTMB2 <- mckeon.glmmTMB1 %>% update(~ . - (1 | BLOCK) + (SYMBIONT | BLOCK),
  control = glmmTMBControl(
    optimizer = optim,
    optArgs = list(method = "SANN")
  )
)
AICc(mckeon.glmmTMB1, mckeon.glmmTMB2)

IN addition to numerous warnings, there is an error about NA/NaN gradient evaluation. This is a somewhat generic error returned by the nlminb optimiser. It occurs when there is an error in either the objective function or gradient function calculations.

Ultimately, this might suggest a poorly specified or over-fit model. We will drop the random intercept/slope.

Model validation

glmer (lme4)

plot_grid

mckeon.glmer1 %>%
  plot_model(type = "diag") %>%
  plot_grid()

Conclusions:

• this is not outrageously bad

Performance model checking

mckeon.glmer1 %>% performance::check_model()

Conclusions:

• it is always difficult to assess pattern in residuals from binary data
• the Q-Q plot seems ok, although the tails are not ideal
• the homogeneity of variance plot seems to suggest a curved pattern, but then again, this is binary data. DHARMa residuals should be easier to diagnose.

DHARMa residuals

mckeon.resid <- mckeon.glmer1 %>% simulateResiduals(plot = TRUE)

mckeon.resid %>% testZeroInflation()

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0428, p-value = 0.936
## alternative hypothesis: two.sided

Conclusions:

• all diagnostics and tests suggest that the various assumptions have been satisified.

glmmTMB (glmmTMB)

plot_grid

mckeon.glmmTMB1 %>%
  plot_model(type = "diag") %>%
  plot_grid()

Conclusions:

• this is not outrageously bad

Performance model checking

mckeon.glmmTMB1 %>% performance::check_model()

Conclusions:

• it is always difficult to assess pattern in residuals from binary data
• the Q-Q plot seems ok, although the tails are not ideal
• the homogeneity of variance plot seems to suggest a curved pattern, but then again, this is binary data. DHARMa residuals should be easier to diagnose.

DHARMa residuals

mckeon.resid <- mckeon.glmmTMB1 %>% simulateResiduals(plot = TRUE)

mckeon.resid %>% testZeroInflation()

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.96105, p-value = 0.968
## alternative hypothesis: two.sided

mckeon.resid %>% testDispersion()

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.0036, p-value = 0.928
## alternative hypothesis: two.sided

Conclusions:

• all diagnostics and tests suggest that the various assumptions have been satisfied.

Partial plots

glmer (lme4)

plot_model

mckeon.glmer1 %>% plot_model(type = "eff", show.data = TRUE, jitter = c(0.05, 0))

## $SYMBIONT

allEffects

mckeon.glmer1 %>%
  allEffects() %>%
  plot()

ggpredict

mckeon.glmer1 %>%
  ggpredict() %>%
  plot(add.data = TRUE, jitter = c(0.1, 0))

## $SYMBIONT

ggemmeans

mckeon.glmer1 %>%
  ggemmeans(~SYMBIONT) %>%
  plot(add.data = TRUE, jitter = c(0.1, 0))

glmmTMB (glmmTMB)

plot_model

mckeon.glmmTMB1 %>% plot_model(type = "eff", show.data = TRUE, jitter = c(0.05, 0))

## $SYMBIONT

allEffects

mckeon.glmmTMB1 %>%
  allEffects() %>%
  plot()

ggpredict

mckeon.glmmTMB1 %>%
  ggpredict() %>%
  plot(add.data = TRUE, jitter = c(0.05, 0))

## $SYMBIONT

ggemmeans

mckeon.glmmTMB1 %>%
  ggemmeans(~SYMBIONT) %>%
  plot(add.data = TRUE, jitter = c(0.05, 0))

Model investigation / hypothesis testing

glmer (lme4)

summary

mckeon.glmer1 %>% summary()

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: PREDATION ~ SYMBIONT + (1 | BLOCK)
##    Data: mckeon
## 
##      AIC      BIC   logLik deviance df.resid 
##     70.7     82.6    -30.4     60.7       75 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -23.2076  -0.1730   0.0976   0.2980   0.8944 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  BLOCK  (Intercept) 11.81    3.437   
## Number of obs: 80, groups:  BLOCK, 10
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)   
## (Intercept)       5.096      1.811   2.814  0.00490 **
## SYMBIONTcrabs    -3.842      1.465  -2.623  0.00871 **
## SYMBIONTshrimp   -4.431      1.551  -2.856  0.00429 **
## SYMBIONTboth     -5.599      1.724  -3.247  0.00116 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) SYMBIONTc SYMBIONTs
## SYMBIONTcrb -0.625                    
## SYMBIONTshr -0.652  0.736             
## SYMBIONTbth -0.687  0.732     0.767

Conclusions:

• the coefficients are presented on a logit scale. Whilst this is not relavant for the purpose of inference testing, it does make it difficult to interpret the coefficients.
• if we exponentiate the coefficients (\(log(\frac{\rho}{1-\rho})\) -> \(\frac{\rho}{1-\rho}\)), they will be presented on a odds ratio scale, and thus:
  • the intercept (none symbionts) will be 163.33. That is, corals without a symbiont are 163.33 times more likely to be preditated on than not predated on. The odds of predation in this the absence of symbionts is 163.33:1.
  • in the presence of a crab symbiont, the odds of being predated on are only 0.02 times that of the none symbiont group. That is, in the presence of a crab symbiont, the odds of predation decline by 98%.
  • in the presence of a shrimp symbiont, the odds of being predated on are only 0.01 times that of the none symbiont group. That is, in the presence of a shrimp symbiont, the odds of predation decline by 99%.
  • in the presence of both crab and shrimp symbionts, the odds of being predated on are only 0 times that of the none symbiont group. That is, in the presence of both crab and shrimp symbiont, the odds of predation decline by 100%.
• if we backtransform the intercept full to the response scale (probability scale), (\(log(\frac{\rho}{1-\rho})\) -> \(\rho\)), the intercept is interpreted as the probability that corals will be predated in the absence of of symbionts is 1

tidy

mckeon.glmer1 %>% tidy(effect = "fixed", conf.int = TRUE, infer = c(TRUE, TRUE), exponentiate = TRUE)

mckeon.glmer1 %>%
  tidy(effect = "fixed", conf.int = TRUE, infer = c(TRUE, TRUE), exponentiate = TRUE) %>%
  kable()

effect	term	estimate	std.error	statistic	p.value	conf.low	conf.high
fixed	(Intercept)	163.3256181	295.7986843	2.813623	0.0048987	4.6929398	5684.1252579
fixed	SYMBIONTcrabs	0.0214459	0.0314110	-2.623280	0.0087088	0.0012151	0.3785015
fixed	SYMBIONTshrimp	0.0119023	0.0184657	-2.856063	0.0042893	0.0005689	0.2490134
fixed	SYMBIONTboth	0.0037016	0.0063822	-3.247315	0.0011650	0.0001261	0.1086478

Conclusions:

• the coefficients are presented on a logit scale. Whilst this is not relevant for the purpose of inference testing, it does make it difficult to interpret the coefficients.
• if we exponentiate the coefficients (\(log(\frac{\rho}{1-\rho})\) -> \(\frac{\rho}{1-\rho}\)), they will be presented on a odds ratio scale, and thus:
  • the intercept (none symbionts) will be 163.33. That is, corals without a symbiont are 163.33 times more likely to be predated on than not predated on. The odds of predation in this the absence of symbionts is 163.33:1.
  • in the presence of a crab symbiont, the odds of being predated on are only 0.02 times that of the none symbiont group. That is, in the presence of a crab symbiont, the odds of predation decline by 98%.
  • in the presence of a shrimp symbiont, the odds of being predated on are only 0.01 times that of the none symbiont group. That is, in the presence of a shrimp symbiont, the odds of predation decline by 99%.
  • in the presence of both crab and shrimp symbionts, the odds of being predated on are only 0 times that of the none symbiont group. That is, in the presence of both crab and shrimp symbiont, the odds of predation decline by 100%.
• if we back-transform the intercept full to the response scale (probability scale), (\(log(\frac{\rho}{1-\rho})\) -> \(\rho\)), the intercept is interpreted as the probability that corals will be predated in the absence of of symbionts is 1

tab_model

# warning this is only appropriate for html output
mckeon.glmer1 %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)

	PREDATION
Predictors	Odds Ratios	std. Error	CI	p
(Intercept)	163.33	295.80	4.69 – 5684.13	0.005
SYMBIONT [crabs]	0.02	0.03	0.00 – 0.38	0.009
SYMBIONT [shrimp]	0.01	0.02	0.00 – 0.25	0.004
SYMBIONT [both]	0.00	0.01	0.00 – 0.11	0.001
Random Effects
σ2	3.29
τ00 BLOCK	11.81
ICC	0.78
N BLOCK	10
Observations	80
Marginal R2 / Conditional R2	0.228 / 0.832
AIC	70.706

glmmTMB (glmmTMB)

summary

mckeon.glmmTMB1 %>% summary()

##  Family: binomial  ( logit )
## Formula:          PREDATION ~ SYMBIONT + (1 | BLOCK)
## Data: mckeon
## 
##      AIC      BIC   logLik deviance df.resid 
##     62.9     74.8    -26.4     52.9       79 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  BLOCK  (Intercept) 17.83    4.223   
## Number of obs: 80, groups:  BLOCK, 10
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       4.420      1.782   2.481 0.013109 *  
## SYMBIONTcrabs    -3.317      1.313  -2.526 0.011534 *  
## SYMBIONTshrimp   -3.956      1.416  -2.793 0.005218 ** 
## SYMBIONTboth     -5.152      1.565  -3.291 0.000997 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conclusions:

• the coefficients are presented on a logit scale. Whilst this is not relevant for the purpose of inference testing, it does make it difficult to interpret the coefficients.
• if we exponentiate the coefficients (\(log(\frac{\rho}{1-\rho})\) -> \(\frac{\rho}{1-\rho}\)), they will be presented on a odds ratio scale, and thus:
  • the intercept (none symbionts) will be 83.1. That is, corals without a symbiont are 83.1 times more likely to be predated on than not predated on. The odds of predation in this the absence of symbionts is 83.1:1.
  • in the presence of a crab symbiont, the odds of being predated on are only 0.04 times that of the none symbiont group. That is, in the presence of a crab symbiont, the odds of predation decline by 96%.
  • in the presence of a shrimp symbiont, the odds of being predated on are only 0.02 times that of the none symbiont group. That is, in the presence of a shrimp symbiont, the odds of predation decline by 98%.
  • in the presence of both crab and shrimp symbionts, the odds of being predated on are only 0.01 times that of the none symbiont group. That is, in the presence of both crab and shrimp symbiont, the odds of predation decline by 99%.
• if we back-transform the intercept full to the response scale (probability scale), (\(log(\frac{\rho}{1-\rho})\) -> \(\rho\)), the intercept is interpreted as the probability that corals will be predated in the absence of of symbionts is 0.99

tidy

mckeon.glmmTMB1 %>% tidy(effect = "fixed", conf.int = TRUE, infer = c(TRUE, TRUE))

mckeon.glmmTMB1 %>% tidy(effect = "fixed", conf.int = TRUE, infer = c(TRUE, TRUE), exponentiate = TRUE)

Conclusions:

• the coefficients are presented on a logit scale. Whilst this is not relevant for the purpose of inference testing, it does make it difficult to interpret the coefficients.
• if we exponentiate the coefficients (\(log(\frac{\rho}{1-\rho})\) -> \(\frac{\rho}{1-\rho}\)), they will be presented on a odds ratio scale, and thus:
  • the intercept (none symbionts) will be 83.1. That is, corals without a symbiont are 83.1 times more likely to be predated on than not predated on. The odds of predation in this the absence of symbionts is 83.1:1.
  • in the presence of a crab symbiont, the odds of being predated on are only 0.04 times that of the none symbiont group. That is, in the presence of a crab symbiont, the odds of predation decline by 96%.
  • in the presence of a shrimp symbiont, the odds of being predated on are only 0.02 times that of the none symbiont group. That is, in the presence of a shrimp symbiont, the odds of predation decline by 98%.
  • in the presence of both crab and shrimp symbionts, the odds of being predated on are only 0.01 times that of the none symbiont group. That is, in the presence of both crab and shrimp symbiont, the odds of predation decline by 99%.
• if we back-transform the intercept full to the response scale (probability scale), (\(log(\frac{\rho}{1-\rho})\) -> \(\rho\)), the intercept is interpreted as the probability that corals will be predated in the absence of of symbionts is 1

tab_model

# warning this is only appropriate for html output
mckeon.glmmTMB1 %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)

	PREDATION
Predictors	Odds Ratios	std. Error	CI	p
(Intercept)	83.10	148.05	2.53 – 2730.03	0.013
SYMBIONT [crabs]	0.04	0.05	0.00 – 0.48	0.012
SYMBIONT [shrimp]	0.02	0.03	0.00 – 0.31	0.005
SYMBIONT [both]	0.01	0.01	0.00 – 0.12	0.001
Random Effects
σ2	3.29
τ00 BLOCK	17.83
ICC	0.84
N BLOCK	10
Observations	80
Marginal R2 / Conditional R2	0.149 / 0.867
AIC	62.861

Further analyses

In addition to comparing each of the symbiont types against the control group of no symbionts, it might be interesting to investigate whether there are any differences between the predation protection provided by crabs and shrimp, as well as whether having both crabs and shrimp symbionts is different to only a single symbiont type.

These contrasts can be explored via specific contrasts.

SYMBIONT	Crab vs Shrimp	One vs Both	None vs Symbiont
none	0	0	1
crab	1	1/2	-1/3
shrimp	-1	1/2	-1/3
both	0	-1	-1/3

glmer (lme4)

planned contrasts

cmat <- cbind(
  crab_vs_shrimp = c(0, 1, -1, 0),
  one_vs_both = c(0, 1 / 2, 1 / 2, -1),
  symbiont = c(1, -1 / 3, -1 / 3, -1 / 3)
)
round(crossprod(cmat), 1)

##                crab_vs_shrimp one_vs_both symbiont
## crab_vs_shrimp              2         0.0      0.0
## one_vs_both                 0         1.5      0.0
## symbiont                    0         0.0      1.3

## all contrasts orthogonal
mckeon.glmer1 %>%
  emmeans(~SYMBIONT, type = "response") %>%
  contrast(method = list(SYMBIONT = cmat)) %>%
  summary(infer = TRUE)

## Or on an absolute difference scale
mckeon.glmer1 %>%
  emmeans(~SYMBIONT, link = "link") %>%
  regrid() %>%
  contrast(method = list(SYMBIONT = cmat)) %>%
  summary(infer = TRUE)

Conclusions:

• there is no inferential difference in predation risk between crab and shrimp symbionts
• there is no inferential difference in predation risk between one symbiont type and both symbionts.
• there is a significant reduction in predation risk associated with having at least one symbiont compared to having no symbionts.

\(R^2\)

mckeon.glmer1 %>% r.squaredGLMM()

##                   R2m       R2c
## theoretical 0.2282003 0.8318487
## delta       0.2166278 0.7896639

## The delta mehod can be used with for all distributions and link
## functions, jwhile lognormal approximation and trigamma function are
## limited to distributions with logarithmic link. Trigamma-estimate
## is recommended whenever available. Additionally, for binomial
## distributions, theoretical variances exist specific for each link
## function distribution.
mckeon.glmer1 %>% performance::r2_nakagawa()

## # R2 for Mixed Models
## 
##   Conditional R2: 0.832
##      Marginal R2: 0.228

glmmTMB (glmmTMB)

planned contrasts

cmat <- cbind(
  crab_vs_shrimp = c(0, 1, -1, 0),
  one_vs_both = c(0, 1 / 2, 1 / 2, -1),
  symbiont = c(1, -1 / 3, -1 / 3, -1 / 3)
)
round(crossprod(cmat), 1)

##                crab_vs_shrimp one_vs_both symbiont
## crab_vs_shrimp              2         0.0      0.0
## one_vs_both                 0         1.5      0.0
## symbiont                    0         0.0      1.3

## all contrasts orthogonal
## On the link scale
mckeon.glmmTMB1 %>%
  emmeans(~SYMBIONT, type = "link") %>%
  contrast(method = list(SYMBIONT = cmat)) %>%
  summary(infer = TRUE)

## On the ratio (fold) scale
mckeon.glmmTMB1 %>%
  emmeans(~SYMBIONT, type = "response") %>%
  contrast(method = list(SYMBIONT = cmat)) %>%
  summary(infer = TRUE)

## On an absolute difference scale
mckeon.glmmTMB1 %>%
  emmeans(~SYMBIONT, type = "link") %>%
  regrid() %>%
  contrast(method = list(SYMBIONT = cmat)) %>%
  summary(infer = TRUE)

Conclusions:

• there is no inferential difference in predation risk between crab and shrimp symbionts
• there is no inferential difference in predation risk between one symbiont type and both symbionts.
• there is a significant reduction in predation risk associated with having at least one symbiont compared to having no symbionts.

\(R^2\)

mckeon.glmmTMB1 %>% r.squaredGLMM()

##                   R2m       R2c
## theoretical 0.1489321 0.8674549
## delta       0.1437647 0.8373578

## The delta method can be used with for all distributions and link
## functions, while lognormal approximation and trigamma function are
## limited to distributions with logarithmic link. Trigamma-estimate
## is recommended whenever available. Additionally, for binomial
## distributions, theoretical variances exist specific for each link
## function distribution.
mckeon.glmmTMB1 %>% performance::r2_nakagawa()

## # R2 for Mixed Models
## 
##   Conditional R2: 0.867
##      Marginal R2: 0.149

Summary figure

glmer (lme4)

mckeon.glmer1 %>%
  emmeans(~SYMBIONT, type = "response") %>%
  as.data.frame() %>%
  ggplot(aes(y = prob, x = SYMBIONT)) +
  geom_pointrange(aes(ymin = asymp.LCL, ymax = asymp.UCL))

glmmTMB (glmmTMB)

mckeon.glmmTMB1 %>%
  emmeans(~SYMBIONT, type = "response") %>%
  as.data.frame() %>%
  ggplot(aes(y = prob, x = SYMBIONT)) +
  geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL))

References

Mckeon, Seabird, Adrian Stier, Shelby Mcilroy, and Benjamin Bolker. 2012. “Multiple Defender Effects: Synergistic Coral Defense by Mutualist Crustaceans.” Oecologia 169 (February): 1095–1103. https://doi.org/10.1007/s00442-012-2275-2.